/*
 * TRISpH.cpp
 *
 *  Created on: 16 Aug 2011
 *      Author: Allan
 */

#include "TRISpH.h"

// C++ includes
#include <iomanip>
#include <iostream>
#include <memory>
#include <string>
#include <vector>
using namespace std;

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include "../../Activity/AqueousActivity.h"
#include "../../Applications/Equilibrator.h"
#include "../../ElectrolyteSolution.h"
#include "../../Multiphase.h"
#include "../../Reaction.h"
#include "../../ReactionSystem.h"
#include "../../Utilities/PhaseUtils.h"
#include "../../Utilities/Range.h"
#include "../../Utilities/ReactionUtils.h"
#include "PitzerTris.h"

const Multiphase
CreateTRISMultiphase();

const ReactionSystem
CreateTRISReactionSystem(const SpeciesMap& speciesMap);

const Equilibrator
CreateTRISEquilibrator();

void TRISpH(double T, double P)
{
	Equilibrator equilibrator = CreateTRISEquilibrator();
	
	cout << setw(15) << left << "mTRIS";
	cout << setw(15) << left << "mNaCl";
	cout << setw(15) << left << "pH";
	cout << setw(15) << left << "ln(gTRIS)";
	cout << endl;

	BOOST_FOREACH(double mTRIS, Range(0.0, 0.02, 0.005)) 
	{
		BOOST_FOREACH(double mNaCl, Range(0.0, 5.0, 1.0))
		{
			equilibrator.Reset();
			
			equilibrator.AddCondition("ChargeBalanced");
			equilibrator.AddCondition("FixedMass", "H2O", 55.508);
			equilibrator.AddCondition("FixedTotalMass", "Na[+]", mNaCl);
			equilibrator.AddCondition("FixedTotalMass", "Cl[-]", mNaCl);
			equilibrator.AddCondition("FixedTotalMass", "TRIS", mTRIS);
			
			SystemState state = equilibrator.Equilibrate(T, P);
			
			cout << setw(15) << left << mTRIS;
			cout << setw(15) << left << mNaCl;
			cout << setw(15) << left << state.Get("pH");
			cout << setw(15) << left << std::log(state.Get("ActivityCoefficient", "TRIS"));
			cout << endl;
		}
		
		cout << endl;
	}
}

void TRISpHTRange(const vector<double>& rangeT, double P, double mNaCl, double mTRIS)
{
	Equilibrator equilibrator = CreateTRISEquilibrator();
	
	cout << setw(15) << left << "T";
	cout << setw(15) << left << "pH";
	cout << setw(15) << left << "ln(gTRIS)";
	cout << endl;

	cout << fixed << setprecision(5);
	
	BOOST_FOREACH(double T, rangeT) 
	{
		equilibrator.Reset();

		equilibrator.AddCondition("ChargeBalanced");
		equilibrator.AddCondition("FixedMass", "H2O", 55.508);
		equilibrator.AddCondition("FixedTotalMass", "Na[+]", mNaCl);
		equilibrator.AddCondition("FixedTotalMass", "Cl[-]", mNaCl);
		equilibrator.AddCondition("FixedTotalMass", "TRIS", mTRIS);

		SystemState state = equilibrator.Equilibrate(T, P);

		cout << setw(15) << left << T;
		cout << setw(15) << left << state.Get("pH");
		cout << setw(15) << left << std::log(state.Get("ActivityCoefficient", "TRIS"));
		cout << endl;
	}
}

void NoTRISpHTRange(const vector<double>& rangeT, double P, double mNaCl)
{
	System system({"H2O", "NaCl(aq)"});
	
	Equilibrator equilibrator(system);
	
	cout << setw(20) << left << "T";
	cout << setw(20) << left << "pH";
	cout << setw(20) << left << "log(K):OH[-]";
	cout << setw(20) << left << "log(K):NaCl(aq)";
	cout << setw(20) << left << "log(K):HCl(aq)";
	cout << setw(20) << left << "log(K):NaOH(aq)";
	cout << endl;

	BOOST_FOREACH(double T, rangeT) 
	{
		equilibrator.Reset();

		equilibrator.AddCondition("ChargeBalanced");
		equilibrator.AddCondition("FixedMass", "H2O", 55.508);
		equilibrator.AddCondition("FixedTotalMass", "Na[+]", mNaCl);
		equilibrator.AddCondition("FixedTotalMass", "Cl[-]", mNaCl);

		SystemState state = equilibrator.Equilibrate(T, P);

		cout << setw(20) << left << T;
		cout << setw(20) << left << state.Get("pH");
		cout << setw(20) << left << std::log10(state.Get("K", "OH[-]"));
		cout << setw(20) << left << std::log10(state.Get("K", "NaCl(aq)"));
		cout << setw(20) << left << std::log10(state.Get("K", "HCl(aq)"));
		cout << setw(20) << left << std::log10(state.Get("K", "NaOH(aq)"));
		cout << endl;
	}
}

const Multiphase CreateTRISMultiphase()
{
	vector<string> aqueousSpecies = {"H2O", "Na[+]", "Cl[-]", "H[+]", "TRIS", "OH[-]", "NaCl(aq)", "HCl(aq)", "NaOH(aq)", "TRISH[+]"};
	
	ElectrolyteSolution sol(aqueousSpecies);
	
	AqueousActivity aqueousActivity(aqueousSpecies);
	
	aqueousActivity.SetActivityCoefficient("H[+]", PitzerHydronModel(sol));
	aqueousActivity.SetActivityCoefficient("TRIS", PitzerTRISModel(sol));
	aqueousActivity.SetActivityCoefficient("TRISH[+]", PitzerTRISHModel(sol));
	
	Multiphase multiphase({CreateAqueousPhase(aqueousSpecies, aqueousActivity)});
	
	return multiphase;
}

const ReactionSystem CreateTRISReactionSystem(const SpeciesMap& speciesMap)
{
	ReactionEquation equationTRISH(speciesMap);
	
	equationTRISH.Set("TRISH[+]", {"TRIS", "H[+]"}, {1.0, 1.0});
	
	EquilibriumConstant K(EqulibriumConstantTRISH);
	
	Reaction reactionTRISH(equationTRISH, K);
	
	vector<Reaction> reactions;
	
	reactions.push_back(CreateReaction("OH[-]", speciesMap));
	reactions.push_back(CreateReaction("NaCl(aq)", speciesMap));
	reactions.push_back(CreateReaction("HCl(aq)", speciesMap));
	reactions.push_back(CreateReaction("NaOH(aq)", speciesMap));
	reactions.push_back(reactionTRISH);
	
	ReactionSystem eReactions(reactions);
	
	return eReactions;
}

const Equilibrator CreateTRISEquilibrator()
{
	Multiphase multiphase = CreateTRISMultiphase();
	
	ReactionSystem eReactions = CreateTRISReactionSystem(multiphase.GetSpeciesMap());
	
	System system(multiphase, eReactions);
	
	Equilibrator equilibrator(system);
	
	return equilibrator;
}
